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Abstract 

Starting from the ordered state, we investigate the short-time behaviour of the 
hard-disk model. For the positional order, we determine the critical exponents rj 
and z from the dynamic relaxation of the order parameter and the cumulant with 
molecular dynamics simulations. The results are compared with previous Monte 
Carlo (MC) simulations. The bond orientational order is studied with MC dynamics. 
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1 Introduction 



For a long time it was believed that universal scaling behaviour can be found 
only in the long-time regime. Therefore, numerical simulations were performed 
in the thermodynamic equilibrium. Such simulations in vicinity of the critical 
point are affected by the critical slowing down. However, recently Janssen, 
Schaub and Schmittmann [1] showed that universality exists already in the 
early time of the evolution. They discovered that a system with non-conserved 
order parameter and energy (model A) quenched from a high temperature 
state to the critical temperature shows universal short-time behaviour already 
after a microscopic time scale tmic- Starting from an unordered state with 
a small value of the order parameter mo, the order increases with a power 
law M{t) ~ rriQt^, where ^ is a new dynamic exponent. A number of Monte 
Carlo (MC) investigations [2-4] support this short-time behaviour. These sim- 
ulations can be also used to calculate the conventional (static and dynamic) 
exponents as well as the critical point [5]. This may eliminate critical slowing 
down, since the simulations are performed in the short-time regime. 

First simulations of the dynamic relaxation in the short-time regime started 
from an unordered state. However, short-time dynamical scaling can be also 
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found starting from the ordered state (M(t = 0) = 1). There exist no ana- 
lytical calculations for this situation, but several MC simulations were done 
[6,7,2]. Also, all critical exponents, except for the new exponent 6', can be cal- 
culated starting from the ordered state. Up to now, simulations of the dynamic 
relaxation have only been performed with MC dynamics. Normally, molecular 
dynamics (MD) simulations can cause ergodicy problems, since the energy of 
the system is conserved. However, for the two-dimensional hard-disk model 
the potential energy of the allowed configurations does not depend on the 
positions of the particles, but is constant. Therefore, the restriction of energy 
conservation does not lead to a reduction of possible configurations. 

The nature of the two-dimensional melting transition is a longstanding puzzle 
[8,9]. The Kosterlitz-Thouless-Halperin-Nelson- Young (KTHNY) theory [10] 
predicts two continuous phase transitions. The first transition (dislocation un- 
binding) at the melting temperature Tm transforms the solid with quasi-long- 
range positional order and long-range orientational order into a new hexatic 
phase that posses short-range positional order and quasi-long-range orienta- 
tional order. The disclination unbinding transition at T; transforms this hex- 
atic phase into an isotropic phase in which the positional and orientational 
order are short-range. There are several other theoretical scenarios for the 
melting transition in two dimensions [8]. Most of them predict a first-order 
phase transition from the solid to the isotropic phase with a coexistence region 
instead of a hexatic phase. 

Even for the simple hard-disk system no consensus about the nature of the 
melting transition has been established. A large number of simulations of the 
two-dimensional hard-disk model in the thermodynamic equilibrium have been 
performed. A melting transition was first seen in a computer study by Alder 
and Wainwright [11]. They investigated 870 disks with MD methods (constant 
number of particles N, volume V and energy E) and found the transition to 
be first order. However, the results of such a small system are affected by large 
finite-size effects. Recent simulations used MC techniques either in the NVT 
ensemble (constant volume) [12-15] or the NpT ensemble (constant pressure) 
[16,17]. Unfortunately, the results of these simulations are not compatible. 

In this article, we study the short-time behaviour of the two-dimensional hard- 
disk model starting from the ordered state (perfect crystal). First we examine 
the positional order parameter ippos{t) and the cumulant Upos{t) with MD 
simulations. The power law behaviour of these observables is used to deter- 
mine the critical exponents 77 and z. The results are compared with those of 
a previous MC simulation [18]. In the second part we study the bond orien- 
tational order parameter ib^lt) and the cumulant UQ{t) with MC dynamics. 
All simulations are performed in a rectangular box with ratio 2 : \/3, which is 
necessary for the ordered state, and with periodic boundary conditions. The 
disk diameter is set equal to one. 
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2 Positional order 



The positional order parameter ■0pos can be computed via 



V'pos(0 



N 



-^exp(iGTl(i)) 

-'^ i=l 



(1) 



where G denotes a reciprocal lattice vector and fi{t) is the position of particle 
i at time t. G has a magnitude of 27r/a, where a = ^2/(\/3p) is the average 
lattice spacing. The direction of G is fixed to that of a reciprocal lattice vector 
of the perfect crystal (which are unique due to the boundary condition of a 
rectangular box of ratio 2 : y/S). The reason for fixing G is that large crystal 
tilting is not possible since we simulate only the short-time behaviour of the 
system. 

MC simulations of the dynamic relaxation of systems with quasi-long-range 
order were performed for the 6-state clock model [19], the XY model [20,21], 
the fully frustrated XY model [3,21], the quantum XY model [22] and the 
hard-disk model [18]. However, no investigations exist for the relaxation with 
molecular dynamics. Independent of the dynamics, the scaling form of the 
second moment of the order parameter at or above pm is 

^^J{t, L) = b-^^^J{b-H, b-'L) , (2) 



where b denotes the rescahng factor. This leads for sufficiently large L to a 
power law time dependence of the form 

V'pos'(t) ~ r"/^ . (3) 



From a finite-size scaling (FSS) analysis of the time-dependent cumulant 

U,os{t) = 7^^^ - 1 (4) 

one obtains 

tj^t) ^ t'/^ , (5) 

where d is the dimension of the system. One can use ?7pos(^) to determine the 
dynamic critical exponent z and then, with z in hand, the static exponent r] 
from the behaviour of V'pos^l^)- 
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In the following wc simulate the hard-disk model with molecular dynamics at 
the melting density Pm ~ 0.933 [18] (the density is given as usual in reduced 
units) and in the solid phase (p = 1.0) to investigate the time evolution of 
^pos^ and Upos- Since the positional order of the system is quasi-long-range, we 
expect to find a power law behaviour for p > Pm as in the case of the MC study. 
Starting from the ordered state (^pos = 1), i-e. from a perfect crystal with 
lattice spacing a, we release the system to evolve with molecular dynamics. 
Solving exactly the simultaneous classical equations of motion for hard-disks 
means in practice to calculate the next collision point of two particles and their 
new momenta. The initial components of the momenta are chosen randomly 
with a distribution proportional exp(p^/2). Obviously, in the case of a hard- 
core potential a global change in the kinetic energy just leads to a rescaling 
of time. We use systems of 16^, 32^, 64^ and 128^ hard-disks and measure 
the observables up to i = 6 in time intervals of At — 0.04. The number of 
independent data sets ranges from 0(150) for N = 128^ to (9(30 000) for 
N ^ 16^ Statistical errors are calculated by dividing the data into different 
subsamples. Systematic errors are estimated by least square fits of different 
system sizes and different time intervals, i.e. the exponent c of the power law 
behaviour is calculated for each system in small intervals [ti,tj]. For large 
enough systems, c is independent of the system size between tmm and tmax- 
The microscopic time scale tmin depends on the microscopic details, but is 
independent from the system size. In contrast to this, tmax — the time when the 
system starts to show finite-size effects — scales with the number of particles. 

In Fig. 1 we plot the time evolution of ■0pos^ at the melting density pm for 
different system sizes in a double logarithmic scale. The figure shows that the 
power law behaviour starts after a microscopic time scale tmic of approximately 
0.16. For times up to 0.5 the difference between the systems with 64^ and 
128'^ hard-disks is negligible. We use these two systems and time intervals 
of i = [0.2 ... 0.8] or smaller for the determination of the critical exponents. 
Power law fits for the different time intervals and system sizes lead to r)/z — 
0.201(4). The situation in the solid phase at p = 1.0 is very similar and we get 
rj/z = 0.0695(25). To determine z independently, we also measure the time 
evolution of the cumulant t/pos- As before the behaviour can be well described 
by a power law ansatz. From the slope we get z — 1.04(3) at pm, while the 
analysis for p = 1.0 yields z = 1.06(3). Thus we get r] = 0.194(9) at p = 0.933 
and f] = 0.0737(49) at p = 1.0, respectively. The values of rj coincide with the 
results from the dynamic relaxation with MC dynamics [18] (using the same 
methods) within statistical errors, as can be seen from Table 1. For comparison 
we also show the results obtained from conventional FSS [18]. The value of the 
dynamic critical exponent z changes from z ^ 2 for MC dynamics to ^ ~ 1 for 
MD simulations. These are the usual values for local MC updating schemes — 
which can be understood as due to the 'diffusion' of the changes induced by 
the local updating — and MD simulations, respectively. The time when the 
system starts to shows finite-size effects scales approximately with for MD 
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Fig. 1. Second moment of the positional order parameter V'pos^ as a function of 
time starting from the ordered state at p = 0.933. tppos^{t) was calculated in time 
steps of Ai = 0.04, where we used MD simulations. 

simulations (as can be estimated from Fig. 1) as well as for MC investigations. 

Our results show that the determination of the critical exponents in the short- 
time region can be done with molecular dynamics at least in the case of hard- 
core potentials. The disadvantage of the MD simulations compared to the 
MC investigations with the Metropolis algorithm is the scaling of CPU time 
with the number of particles. In the Metropolis algorithm the CPU time for a 
sweep is proportional to the number of particles N, i.e. a single step of a disk 
is independent of A^. In the MD simulations the number of collisions during a 
fixed time interval At scales also with N. However, the search for the next two 

Table 1 



The critical exponents z and rj determined from the short-time behaviour of C/pos(i) 
and V'pos^(^) with MD and MC dynamics [16] and the value of rj measured with FSS 
methods [16]. 





MD short-time 


MC short-time 


FSS 


p 


z rj 


z rj 


V 


0.933 


1.04(3) 0.194(9) 


2.01(2) 0.199(3) 


0.200(2) 


1.0 


1.06(3) 0.0737(49) 


2.06(4) 0.0794(29) 


0.0791(6) 
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colliding disks is not independent of A^. The CPU time for a naive algorithm 
scales with A^^, while an improved version yields a factor A^. Therefore, also the 
improved MD algorithm leads to a poorer scaling behaviour as the MC version. 
However, the MD simulations not only enlarge the knowledge of dynamic 
relaxation but perhaps also offer a possibility to compare the results with 
experimental investigations. 



3 Bond orientational order 

The orientational order of the hard-disk system can be described by the bond 
orientational order parameter 

^6 = T7 TT 13 exp(6 i . (6) 
i=l ^^l j=l 

The sum on j is over the iVj neighbours of the particle i and 9ij is the angle 
between the particles i and j and an arbitrary but fixed reference axis. Two 
particles are defined as neighbours, if the distance is less than 1.4 times the 
average lattice spacing a. This definition is computationally less expensive 
than the precise determination with the Voronoi construction [23]. 

We examine the bond orientational order of the hard-disk model in the liquid 
regime at p = 0.885, in vicinity of the transition point pi ~ 0.899 [15] at 
p — 0.9 and p — 0.905 and in the solid phase at p = 0.94. We measure the 
second moment of the order parameter i/jq^ and the cumulant Ue ets & function 
of time, where the MC dynamics is given by the Metropolis algorithm [24]. 
The new positions of the particles are chosen with equal probability within a 
circle centered about its original position. As before, we start the relaxation 
process from the perfect ordered crystal {ipe = i^pos = !)• We use systems of 
64^ and 128^ hard-disks and measure up to 10 000 MC sweeps. 

In case of a KTHNY-like scenario pi is the beginning of the hexatic phase, i.e., 
the lower bond of the critical line which ends at pm- Therefore, we expect a 
power law behaviour similar to Eq. (3) between pi and pm assuming a scaling 
behaviour of the form (2) for the positional order parameter. The value of the 
critical exponent 1]^ (at p = pi) is predicted with 1/4 [10] and was measured 
with 0.251(36) [15], while the value of the dynamic critical exponent z for the 
local Metropolis algorithm is normally about two. For a conventional weak 
first-order phase transition the behaviour of the positional order parameter 
should be also approximately power-like. 

Figure 2 shows the time evolution of t/'e^ for the different densities. Obviously, 
the time dependence in vicinity of pi can not be described by a simple power 
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Fig. 2. Time evolution of the second moment of the bond orientational order pa- 
rameter starting from the ordered state at p = 0.885, 0.900, 0.905 and 0.940. 
The MC dynamics is given by the Metropohs algorithm. 

law behaviour as expected. Therefore, the scahng form of the second moment 
of the order parameter T/^g^ is not given by (2). The behaviour of ip%'{t) at 
p = 0.9 and 0.905 is also not consistent with a simple weak first-order phase 
transition. The behaviour in the fluid phase {p = 0.885) is incompatible with 
both scenarios, since it is not of the form t"^^/^ exp(— t/^t) [25]. To rule out 
that the non power law behaviour is just an effect coming from our method 
determing neighbours, we also make simulations using the precise Voronoi 
definition. The result for ipQ^{t) at p = 0.9 is visualized in Fig. 3. For small 
times the difference between both definitions is negligible. At larger times the 
number of disclinations increase and causes errors in case of using the distance 
for the determination of neighbours. Therefore, the deviation between the 
different values of tpQ grows if the time increases. However, also in case of 
using the Voronoi construction we find no power law dependence. 

A possible explanation for the behaviour of ipe^i'^) could be a mixing of the 
order parameters^. In the simplest case we expect that the behaviour should 
be described by the sum of two power law functions. However, we find that 
the the curves at p = 0.9 and p = 0.905 could not be fitted very well by such 
an ansatz. Nevertheless, we use the almost linear behaviour in the intervals 
ti = [10, 100] and t2 = [1000, 10 000] to determine the values of the exponents, 

^ A mixing of order parameters was also found for the Ashkin- Teller model [26]. 
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Fig. 3. Time evolution of the second moment of the bond orientational order pa- 
rameter starting from the ordered state at p = 0.9. The soHd curve is the same 
as in Fig. 2, i.e. neighbours are defined with the distance of the particles. The dashed 
line shows the behaviour that we obtained using the Voronoi construction for the 
determination of neighbours. 



i.e. we examine if the behaviour in one of these intervals is consistent with a 
power law coming from the bond orientational order parameter. At p = 0.9 we 
get from ip^^if) the exponent Ci = 0.124 and 0.0264, respectively. The slope 
of the cumulant UQ{t) (which is shown in Fig. 4) gives the exponent c\j. This 
yields z = 2.12, r) = 0.262 (for ti = [10, 100]) and z = 4.95, r] = 0.130 (for t2 = 
[1000, 10 000]). The first value of r/ is consistent with previous measurements of 
in equilibrium [15], while the value determined in the interval t2 is too small. 
However, also the behaviour oiip^it) in the interval ti is inconsistent with (2) 
since the value is not constant in the solid phase. Thus, neither the behaviour 
in the time interval ti nor in the interval ^2 is compatible with the scaling form 
(2). Additionally, we try to estimate the dynamic critical exponent z from the 
exponential behaviour the the fluid phase. But these measurements give also 
inconsistent results, i.e. we get a value z \. 
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Fig. 4. Time dependent cumulant Uait) at p = 0.9 for N = 64^ and N = 128^ for 
the Metropolis algorithm. 



4 Conclusions 



We investigated the short-time behaviour of the hard-disk model for the po- 
sitional and bond orientational order. The positional order parameter was 
studied with MD simulations at pm and in the solid phase. We determined 
the critical exponent rj as well as the dynamic critical exponent z from the 
power law behaviour of ■0pos^(^) and t/pos(^)- The values of r] are in agreement 
with previous measurements. Our results show that MD simulations can be 
used for the determination of critical exponents from the short-time behaviour. 
The dynamic critical exponent changes from 2; ~ 2 for MC to 2; ~ 1 for MD 
simulations. 



The bond orientational order was studied with conventional MC dynamics. 
The time evolution of the order parameter and the cumulant is not given by a 
simple power law behaviour, i.e. the scaling behaviour is not of the form (2). 
This phenomena could not be explained by a weak first-order phase transition 
or by the way of determining neighbours. A possible mixing of order parameter 
was discussed. 
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